## add 'developer' to packages to be installed from github
x <- c("remotes", "lubridate", "readxl", "pbapply", "viridis", "ggplot2", "kableExtra", "knitr", "formatR", "MASS", "sp", "GGally", "brms", "lme4", "dplyr", "purrr", "forcats", "tidyr", "modelr", "tidybayes", "cowplot", "ggrepel", "posterior", "ggridges", "maRce10/PhenotypeSpace")
aa <- lapply(x, function(y) {
# get pakage name
pkg <- strsplit(y, "/")[[1]]
pkg <- pkg[length(pkg)]
# check if installed, if not then install
if (!pkg %in% installed.packages()[,"Package"]) {
if (grepl("/", y)) remotes::install_github(y, force = TRUE) else
install.packages(y)
}
# load package
a <- try(require(pkg, character.only = T), silent = T)
if (!a) remove.packages(pkg)
})
opts_knit$set(root.dir = "..")
# set evaluation false
opts_chunk$set(fig.width = 10, fig.height = 6, warning = FALSE, message = FALSE, tidy = TRUE)
read_excel_df <- function(...) data.frame(read_excel(...))
# for reading months in english format
sl <- Sys.setlocale(locale = "en_US.UTF-8")
standard_error <- function(x) sd(x)/sqrt(length(x))
cols <- viridis(10, alpha = 0.7)
# print results brms models
summary_brm_model <- function(x, gsub.pattern = NULL, gsub.replacement = NULL, xlab = "Effect size"){
summ <- summary(x)$fixed
fit <- x$fit
betas <- grep("^b_", names(fit@sim$samples[[1]]), value = TRUE)
hdis <- t(sapply(betas, function(y) hdi(fit@sim$samples[[1]][[y]]))
)
summ_table <- data.frame(summ, hdis, iterations = attr(fit, "stan_args")[[1]]$iter, chains = length(attr(fit, "stan_args")))
summ_table <- summ_table[rownames(summ_table) != "Intercept", c("Estimate", "Rhat", "Bulk_ESS", "l.95..CI", "u.95..CI", "iterations", "chains")]
summ_table <- as.data.frame(summ_table)
summ_table$Rhat <- round(summ_table$Rhat, digits = 3)
summ_table$CI_low <- round(unlist(summ_table$l.95..CI), digits = 3)
summ_table$CI_high <- round(unlist(summ_table$u.95..CI), digits = 3)
summ_table$l.95..CI <- summ_table$u.95..CI <- NULL
out <- lapply(betas, function(y) data.frame(variable = y, value = sort(fit@sim$samples[[1]][[y]], decreasing = FALSE)))
posteriors <- do.call(rbind, out)
posteriors <- posteriors[posteriors$variable != "b_Intercept", ]
if (!is.null(gsub.pattern) & !is.null(gsub.replacement))
posteriors$variable <- gsub(pattern = gsub.pattern, replacement = gsub.replacement, posteriors$variable)
gg <- ggplot(data = posteriors, aes(y = variable, x = value, fill = stat(quantile))) +
geom_density_ridges_gradient(quantile_lines = TRUE, quantile_fun = HDInterval::hdi, vline_linetype = 2) +
theme_classic() +
labs(x = xlab, y = "Predictor") +
xlim(range(c(min(summ_table[ , "CI_low"]), max(summ_table[ , "CI_high"])), 0)) +
geom_vline(xintercept = 0, col = "red") +
scale_fill_manual(values = c("transparent", "lightblue", "transparent"), guide = "none") + ggtitle(x$formula)
if (!is.null(gsub.pattern) & !is.null(gsub.replacement))
rownames(summ_table) <- gsub(pattern = gsub.pattern, replacement = gsub.replacement, rownames(summ_table))
summ_table$Rhat <- ifelse(summ_table$Rhat > 1.05, cell_spec(summ_table$Rhat, "html", color ="white", background = "red", bold = TRUE, font_size = 12), cell_spec(summ_table$Rhat, "html"))
signif <- summ_table[,"CI_low"] * summ_table[,"CI_high"] > 0
df1 <- kbl(summ_table, row.names = TRUE, escape = FALSE, format = "html", digits = 3)
df1 <- row_spec(kable_input = df1,row = which(summ_table$CI_low * summ_table$CI_high > 0), background = adjustcolor(cols[9], alpha.f = 0.3))
df1 <- kable_styling(df1, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, font_size = 12)
print(df1)
print(gg)
}
# read ext sel tab calls
sels <- read.csv("./data/processed/tailored_budgie_calls_sel_tab.csv")
# keep only spectrographic parameters
sels <- sels[, c("sound.files", "selec", "duration", "meanfreq", "sd", "freq.median",
"freq.IQR", "time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom",
"mindom", "maxdom", "dfrange", "modindx", "meanpeakf")]
sels$ID <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][1])
sels$month <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][2])
sels$day <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][3])
sels$year <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][4])
sels$date <- paste(sels$day, substr(sels$month, 0, 3), sels$year, sep = "-")
sels$date <- as.Date(sels$date, format = "%d-%b-%Y")
# acoustic measurements
areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week.RDS")
indiv_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_first_week_by_individual.RDS")
indiv_ovlp$treatment <- factor(indiv_ovlp$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
group_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_group_by_week.RDS")
# group_ovlp <- do.call(rbind, lapply(unique(group_ovlp$ID), function(x) { X <-
# group_ovlp[group_ovlp$ID == x, ] X$overlap.to.group <- X$overlap.to.group -
# X$overlap.to.group[X$week == min(as.numeric(as.character(X$week)))]
# X$distance.to.group <- X$distance.to.group - X$distance.to.group[X$week ==
# min(as.numeric(as.character(X$week)))] return(X) }))
group_ovlp$treatment <- factor(group_ovlp$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
df <- as.data.frame(table(sels$ID))
names(df) <- c("ID", "Sample_size")
df <- df[order(df$Sample_size, decreasing = FALSE), ]
kb <- kable(df, row.names = FALSE)
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
print(kb)
| ID | Sample_size |
|---|---|
| 132YGMM | 6 |
| 125YGMM | 12 |
| 160YGHM | 16 |
| 310YGHM | 21 |
| 393YGHM | 25 |
| 303YGHM | 37 |
| 398WBHM | 41 |
| 365YGHM | 44 |
| 399YGLM | 46 |
| 300YGHM | 47 |
| 270YGHM | 64 |
| 407YGHM | 97 |
| 386WBHM | 100 |
| 377WWLM | 108 |
| 367WWNM | 119 |
| 371YYLM | 149 |
| 397YGHM | 175 |
| 378YYLM | 195 |
| 404WBHM | 196 |
| 258YGHM | 223 |
| 389WWLM | 230 |
| 262YGHM | 306 |
| 400YGHM | 360 |
| 388YGLM | 377 |
| 327YYLM | 404 |
| 396YBHM | 455 |
| 403WBHM | 566 |
| 356WBHM | 761 |
| 361YGHM | 767 |
| 402YGHM | 849 |
| 395WBHM | 854 |
| 360YGHM | 900 |
| 390WBHM | 935 |
| 405YBHM | 975 |
| 385YBHM | 1256 |
| 394YBHM | 1693 |
metadat <- read_excel_df("./data/raw/INBREStress_MasterDataSheet_14Nov19.xlsx")
# head(metadat)
sels$ID[sels$ID == "125YGMM"] <- "125YGHM"
sels$ID[sels$ID == "394YBHM"] <- "394WBHM"
# setdiff(sels$ID, metadat$Bird.ID) setdiff(metadat$Bird.ID, sels$ID)
sels$treatment <- sapply(1:nrow(sels), function(x) {
metadat$Treatment[metadat$Bird.ID == sels$ID[x]][1]
})
sels$treatment.room <- sapply(1:nrow(sels), function(x) {
metadat$Treatment.Room[metadat$Bird.ID == sels$ID[x]][1]
})
sels$round <- sapply(1:nrow(sels), function(x) {
metadat$Round[metadat$Bird.ID == sels$ID[x]][1]
})
sels$source.room <- sapply(1:nrow(sels), function(x) {
metadat$Source.Room[metadat$Bird.ID == sels$ID[x]][1]
})
sels$record.group <- sapply(1:nrow(sels), function(x) {
metadat$Record.Group[metadat$Bird.ID == sels$ID[x]][1]
})
# add week
out <- lapply(unique(sels$round), function(x) {
Y <- sels[sels$round == x, ]
min_date <- min(Y$date)
week_limits <- min_date + seq(0, 100, by = 7)
Y$week <- NA
for (i in 2:length(week_limits)) Y$week[Y$date >= week_limits[i - 1] & Y$date <
week_limits[i]] <- i - 1
return(Y)
})
sels <- do.call(rbind, out)
sels$cort.baseline <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
sels$cort.stress <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
sels$stress.response <- sels$cort.stress - sels$cort.baseline
sels$weight <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
sels$breath.count <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.Breath.Count[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.Breath.Count[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
# remove week 5
sels <- sels[sels$week != 5, ]
agg_dat <- aggregate(selec ~ ID + week, data = sels, length)
# compare to week 1 agg_dat$call.count <- sapply(1:nrow(agg_dat), function(x) {
# baseline <- agg_dat$selec[agg_dat$week == 1 & agg_dat$ID == agg_dat$ID[x]] if
# (length(baseline) > 0) change <- agg_dat$selec[x] - baseline else change <-
# agg_dat$selec[x] return(change) } )
# without comparing to week 1
agg_dat$call.count <- sapply(1:nrow(agg_dat), function(x) agg_dat$selec[x])
agg_dat$selec <- NULL
agg_dat$baseline.CORT <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$cort.baseline[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$cort.baseline[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$stress.response <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$stress.response[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$stress.response[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$stress.CORT <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$cort.stress[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$cort.stress[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$weight <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$weight[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$weight[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$breath.rate <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$breath.count[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$breath.count[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$acoustic.area <- sapply(1:nrow(agg_dat), function(x) {
area <- areas_by_week$raref.area[areas_by_week$ID == agg_dat$ID[x] & areas_by_week$week ==
agg_dat$week[x]]
if (length(area) < 1)
area <- NA
return(area)
})
agg_dat$acoustic.distance <- sapply(1:nrow(agg_dat), function(x) {
distance <- indiv_ovlp$distance.to.first.week[indiv_ovlp$ID == agg_dat$ID[x] &
indiv_ovlp$week == agg_dat$week[x]]
if (length(distance) < 1)
distance <- NA
return(distance)
})
agg_dat$acoustic.overlap <- sapply(1:nrow(agg_dat), function(x) {
overlap <- indiv_ovlp$overlap.to.first.week[indiv_ovlp$ID == agg_dat$ID[x] &
indiv_ovlp$week == agg_dat$week[x]]
if (length(overlap) < 1)
overlap <- NA
return(overlap)
})
agg_dat$acoustic.overlap.to.group <- sapply(1:nrow(agg_dat), function(x) {
overlap <- group_ovlp$overlap.to.group[group_ovlp$ID == agg_dat$ID[x] & group_ovlp$week ==
agg_dat$week[x]]
if (length(overlap) < 1)
overlap <- NA
return(overlap)
})
agg_dat$treatment <- sapply(1:nrow(agg_dat), function(x) unique(sels$treatment[sels$ID ==
agg_dat$ID[x]]))
agg_dat$round <- sapply(1:nrow(agg_dat), function(x) unique(sels$round[sels$ID ==
agg_dat$ID[x]]))
breath.count <- stack(metadat[, c("D3.Breath.Count", "D7.Breath.Count", "D14.Breath.Count",
"D21.Breath.Count", "D28.Breath.Count")])
weight <- stack(metadat[, c("D3.Bird.Weight..g.", "D7.Bird.Weight..g.", "D14.Bird.Weight..g.",
"D21.Bird.Weight..g.", "D28.Bird.Weight..g.")])
cort.stress <- stack(metadat[, c("D3.CORT.Stress", "D7.CORT.Stress", "D14.CORT.Stress",
"D21.CORT.Stress", "D28.CORT.Stress")])
cort.baseline <- stack(metadat[, c("D3.CORT.Baseline", "D7.CORT.Baseline", "D14.CORT.Baseline",
"D21.CORT.Baseline", "D28.CORT.Baseline")])
stress <- data.frame(metadat[, c("Bird.ID", "Treatment", "Round", "Treatment.Room")],
week = breath.count$ind, breath.count = breath.count$values, weight = weight$values,
cort.stress = cort.stress$values, cort.baseline = cort.baseline$values, stress.response = cort.stress$values -
cort.baseline$values)
# head(stress)
stress$week <- factor(sapply(strsplit(as.character(stress$week), "\\."), "[[", 1),
levels = c("D3", "D7", "D14", "D21", "D28"))
stress$day <- as.numeric(gsub("D", "", as.character(stress$week)))
stress$round <- paste("Round", stress$Round)
stress$treatment <- factor(stress$Treatment, levels = c("Control", "Medium Stress",
"High Stress"))
# remove 5th week
stress <- stress[stress$week != "D28", ]
stress_l <- lapply(stress$Bird.ID, function(x) {
X <- stress[stress$Bird.ID == x, ]
X$breath.count <- X$breath.count - X$breath.count[X$week == "D3"]
X$weight <- X$weight - X$weight[X$week == "D3"]
X$cort.stress <- X$cort.stress - X$cort.stress[X$week == "D3"]
X$cort.baseline <- X$cort.baseline - X$cort.baseline[X$week == "D3"]
X$stress.response <- X$stress.response - X$stress.response[X$week == "D3"]
return(X)
})
stress <- do.call(rbind, stress_l)
agg_stress <- aggregate(cbind(breath.count, weight, stress.response, cort.baseline) ~
week + treatment, stress, mean)
agg_stress_se <- aggregate(cbind(breath.count, weight, stress.response, cort.baseline) ~
week + treatment, stress, standard_error)
names(agg_stress_se) <- paste(names(agg_stress_se), ".se", sep = "")
agg_stress <- cbind(agg_stress, agg_stress_se[, 3:6])
agg_stress$Week <- 1:4
bs <- 22
gg_breath.count <- ggplot(data = agg_stress, aes(x = Week, y = breath.count, fill = treatment)) +
geom_bar(stat = "identity") + geom_errorbar(aes(ymin = breath.count - breath.count.se,
ymax = breath.count + breath.count.se), width = 0.1) + scale_fill_viridis_d(begin = 0.1,
end = 0.9) + facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Mean change in breath\nrate (breaths/min)",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_weight <- ggplot(data = agg_stress, aes(x = Week, y = weight, fill = treatment)) +
geom_bar(stat = "identity") + geom_errorbar(aes(ymin = weight - weight.se, ymax = weight +
weight.se), width = 0.1) + scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
ncol = 3, scale = "fixed") + labs(y = "Mean change in \nweight (grams)", x = "Week") +
theme_classic(base_size = bs) + theme(legend.position = "none")
gg_cort.baseline <- ggplot(data = agg_stress, aes(x = Week, y = cort.baseline, fill = treatment)) +
geom_bar(stat = "identity") + geom_errorbar(aes(ymin = cort.baseline - cort.baseline.se,
ymax = cort.baseline + cort.baseline.se), width = 0.1) + scale_fill_viridis_d(begin = 0.1,
end = 0.9) + facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Mean change in\nbaseline CORT (ng/mL)",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_stress.response <- ggplot(data = agg_stress, aes(x = Week, y = stress.response,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = stress.response -
stress.response.se, ymax = stress.response + stress.response.se), width = 0.1) +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment, ncol = 3,
scale = "fixed") + labs(y = "Mean change in stress\nresponse CORT (ng/mL)", x = "Week") +
theme_classic(base_size = bs) + theme(legend.position = "none")
cowplot::plot_grid(gg_breath.count, gg_weight, gg_cort.baseline, gg_stress.response)
Models: Predicted physio measure ~ treatment + week (continuous) + IndRandom
Variables (Difference from Week 1): weight, BR, baseline CORT, Stress CORT, Stress Response
responses <- c("baseline.CORT", "stress.response", "stress.CORT", "weight", "breath.rate")
predictors <- c("~ treatment + week + (1|ID) + (1|round)")
formulas <- expand.grid(responses = responses, predictors = predictors, stringsAsFactors = FALSE)
vars_to_scale <- c(responses, "week")
# remove week 1
sub_agg_dat <- agg_dat[agg_dat$week != 1, ]
for (i in vars_to_scale) sub_agg_dat[, vars_to_scale] <- scale(sub_agg_dat[, vars_to_scale])
physio_models <- lapply(1:nrow(formulas), function(x) {
sub_dat <- sub_agg_dat[!is.na(sub_agg_dat[names(sub_agg_dat) == formulas$responses[x]]),
]
sub_dat
mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
iter = 20000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
chains = 4), silent = TRUE)
return(mod)
})
names(physio_models) <- paste(formulas$responses, formulas$predictors)
saveRDS(physio_models, "./data/processed/physiological_response_models.RDS")
physio_models <- readRDS("./data/processed/physiological_response_models.RDS")
ggl <- list()
for (x in 1:length(physio_models)) {
cat(paste("- ", names(physio_models)[x], "\n"))
ggl[[length(ggl) + 1]] <- summary_brm_model(physio_models[[x]], gsub.pattern = "b_treatment|b_",
gsub.replacement = "")
}
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| treatmentHighStress | 0.763 | 1.001 | 10563.66 | 20000 | 4 | 0.114 | 1.417 |
| treatmentMediumStress | 0.180 | 1 | 11354.05 | 20000 | 4 | -0.511 | 0.868 |
| week | 0.000 | 1 | 30967.77 | 20000 | 4 | -0.142 | 0.142 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| treatmentHighStress | -0.493 | 1.001 | 8754.839 | 20000 | 4 | -1.182 | 0.191 |
| treatmentMediumStress | 0.002 | 1.001 | 7335.265 | 20000 | 4 | -0.768 | 0.760 |
| week | -0.123 | 1.001 | 38086.430 | 20000 | 4 | -0.255 | 0.007 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| treatmentHighStress | -0.217 | 1 | 12610.60 | 20000 | 4 | -0.908 | 0.469 |
| treatmentMediumStress | 0.076 | 1 | 13662.12 | 20000 | 4 | -0.697 | 0.837 |
| week | -0.153 | 1 | 37967.80 | 20000 | 4 | -0.285 | -0.020 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| treatmentHighStress | -0.315 | 1 | 10029.64 | 20000 | 4 | -0.966 | 0.343 |
| treatmentMediumStress | -0.126 | 1 | 10793.59 | 20000 | 4 | -0.835 | 0.587 |
| week | -0.344 | 1 | 21762.91 | 20000 | 4 | -0.478 | -0.210 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| treatmentHighStress | 0.145 | 1.001 | 21047.468 | 20000 | 4 | -0.158 | 0.438 |
| treatmentMediumStress | 0.058 | 1.001 | 8702.607 | 20000 | 4 | -0.275 | 0.380 |
| week | -0.785 | 1 | 35693.191 | 20000 | 4 | -0.898 | -0.670 |
names(ggl) <- names(physio_models)
Barplot and effect sizes side-by-side
cowplot::plot_grid(gg_breath.count + theme_classic(base_size = 10) + theme(legend.position = "none"),
ggl[[grep("breath", names(ggl))]])
cowplot::plot_grid(gg_weight + theme_classic(base_size = 10) + theme(legend.position = "none"),
ggl[[grep("weight", names(ggl))]])
cowplot::plot_grid(gg_cort.baseline + theme_classic(base_size = 10) + theme(legend.position = "none"),
ggl[[grep("base", names(ggl))]])
cowplot::plot_grid(gg_stress.response + theme_classic(base_size = 10) + theme(legend.position = "none"),
ggl[[grep("response", names(ggl))]])
Breath rate decreases gradually with time across after the first week
Stress response is higher in “high stress” birds compared to first week
t-SNE
scale_param <- scale(sels[, c("duration", "meanfreq", "sd", "freq.median", "freq.IQR",
"time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom", "mindom",
"maxdom", "dfrange", "modindx", "meanpeakf")])
tsne <- Rtsne(scale_param, dims = 2, perplexity = 30, verbose = FALSE, max_iter = 5000)
saveRDS(tsne, "./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")
tsne <- readRDS("./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")
Y <- as.data.frame(tsne$Y)
names(Y) <- c("TSNE1", "TSNE2")
sels <- data.frame(sels, Y)
sels$treatment <- factor(sels$treatment, levels = c("Control", "Medium Stress", "High Stress"))
ggplot(sels, aes(x = TSNE1, y = TSNE2, col = as.factor(treatment))) + geom_point() +
labs(color = "Treatment") + scale_color_viridis_d(alpha = 0.4) + theme_classic(base_size = 25) +
guides(colour = guide_legend(override.aes = list(size = 10)))
agg_call.count <- aggregate(cbind(call.count, acoustic.overlap.to.group) ~ week +
treatment, agg_dat, mean)
agg_behav <- aggregate(cbind(acoustic.area, acoustic.overlap) ~ week + treatment,
agg_dat, mean)
agg_call.count_se <- aggregate(cbind(call.count, acoustic.overlap.to.group) ~ week +
treatment, agg_dat, standard_error)
agg_behav_se <- aggregate(cbind(acoustic.area, acoustic.overlap) ~ week + treatment,
agg_dat, standard_error)
agg_behav_se <- merge(agg_call.count_se, agg_behav_se, all = TRUE)
names(agg_behav_se) <- paste(names(agg_behav_se), ".se", sep = "")
agg_behav <- merge(agg_call.count, agg_behav, all = TRUE)
agg_behav <- cbind(agg_behav, agg_behav_se[, 3:6])
bs <- 22
agg_behav$treatment <- factor(agg_behav$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
gg_call.count <- ggplot(data = agg_behav, aes(x = week, y = call.count, fill = treatment)) +
geom_bar(stat = "identity") + geom_errorbar(aes(ymin = call.count - call.count.se,
ymax = call.count + call.count.se), width = 0.1) + scale_fill_viridis_d(begin = 0.1,
end = 0.9) + facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Vocal output (number of calls)",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_acoustic.area <- ggplot(data = agg_behav, aes(x = week, y = acoustic.area, fill = treatment)) +
geom_bar(stat = "identity") + geom_errorbar(aes(ymin = acoustic.area - acoustic.area.se,
ymax = acoustic.area + acoustic.area.se), width = 0.1) + scale_fill_viridis_d(begin = 0.1,
end = 0.9) + facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Mean change in \n acoustic area",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_acoustic.overlap <- ggplot(data = agg_behav, aes(x = week, y = acoustic.overlap,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = acoustic.overlap -
acoustic.overlap.se, ymax = acoustic.overlap + acoustic.overlap.se), width = 0.1) +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment, ncol = 3,
scale = "fixed") + labs(y = "Acoustic overlap to self in week 1", x = "Week") +
theme_classic(base_size = bs) + theme(legend.position = "none")
gg_acoustic.overlap.group <- ggplot(data = agg_behav, aes(x = week, y = acoustic.overlap.to.group,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = acoustic.overlap.to.group -
acoustic.overlap.to.group.se, ymax = acoustic.overlap.to.group + acoustic.overlap.to.group.se),
width = 0.1) + scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
ncol = 3, scale = "fixed") + labs(y = "Acoustic overlap to group", x = "Week") +
theme_classic(base_size = bs) + theme(legend.position = "none")
cowplot::plot_grid(gg_call.count, gg_acoustic.area, gg_acoustic.overlap, gg_acoustic.overlap.group)
Model: Predicted behavior ~ treatment + week (continuous) + IndRandom
Variables: # calls, Distance moved from self in first week, Overlap to original acoustic space, Match to group repertoire, Maybe overall size of acoustic space
Do as comparison to week one using rarefacted calls and kernel density
responses <- c("call.count", "acoustic.area", "acoustic.overlap", "acoustic.overlap.to.group")
predictors <- c("~ treatment + week + (1|ID) + (1|round)")
formulas <- expand.grid(responses = responses, predictors = predictors, stringsAsFactors = FALSE)
vars_to_scale <- c(responses, "week")
for (i in vars_to_scale) agg_dat[, vars_to_scale] <- scale(agg_dat[, vars_to_scale])
behav_models <- lapply(1:nrow(formulas), function(x) {
sub_dat <- agg_dat[!is.na(agg_dat[names(agg_dat) == formulas$responses[x]]),
]
# remove week 1
if (!grepl("count|group", formulas$responses[x]))
sub_dat <- sub_dat[sub_dat$week != 1, ]
mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
iter = 50000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
chains = 4), silent = TRUE)
return(mod)
})
names(behav_models) <- paste(formulas$responses, formulas$predictors)
saveRDS(behav_models, "./data/processed/behavioral_response_models.RDS")
behav_models <- readRDS("./data/processed/behavioral_response_models.RDS")
ggl <- list()
for (x in 1:length(behav_models)) {
cat(paste("- ", names(behav_models)[x], "\n \n---\n "))
ggl[[length(ggl) + 1]] <- summary_brm_model(behav_models[[x]], gsub.pattern = "b_treatment|b_",
gsub.replacement = "")
}
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| treatmentHighStress | 0.790 | 1.001 | 11548.09 | 50000 | 4 | 0.213 | 1.380 |
| treatmentMediumStress | 0.398 | 1.001 | 17599.85 | 50000 | 4 | -0.182 | 0.985 |
| week | -0.059 | 1 | 102691.50 | 50000 | 4 | -0.203 | 0.086 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| treatmentHighStress | -0.206 | 1.001 | 10719.216 | 50000 | 4 | -1.059 | 0.624 |
| treatmentMediumStress | -0.708 | 1 | 18329.938 | 50000 | 4 | -1.567 | 0.143 |
| week | 0.035 | 1.003 | 1444.876 | 50000 | 4 | -0.094 | 0.152 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| treatmentHighStress | 0.910 | 1 | 25225.85 | 50000 | 4 | 0.065 | 1.757 |
| treatmentMediumStress | 0.233 | 1 | 25160.95 | 50000 | 4 | -0.699 | 1.151 |
| week | -0.169 | 1 | 76169.60 | 50000 | 4 | -0.307 | -0.032 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| treatmentHighStress | 0.351 | 1.001 | 28237.431 | 50000 | 4 | -0.292 | 1.015 |
| treatmentMediumStress | 0.380 | 1.002 | 3209.631 | 50000 | 4 | -0.280 | 1.046 |
| week | -0.223 | 1.001 | 41111.148 | 50000 | 4 | -0.385 | -0.057 |
names(ggl) <- names(behav_models)
Barplot and effect sizes side-by-side
cowplot::plot_grid(gg_call.count + theme_classic(base_size = 10) + theme(legend.position = "none"),
ggl[[grep("count", names(ggl))]])
cowplot::plot_grid(gg_acoustic.area + theme_classic(base_size = 10) + theme(legend.position = "none"),
ggl[[grep("area", names(ggl))]])
cowplot::plot_grid(gg_acoustic.overlap + theme_classic(base_size = 10) + theme(legend.position = "none"),
ggl[[grep("overlap ~", names(ggl))]])
cowplot::plot_grid(gg_acoustic.overlap.group + theme_classic(base_size = 10) + theme(legend.position = "none"),
ggl[[grep("group", names(ggl))]])
Higher vocal output in “high stress” birds compared to control
Higher acoustic overlap to themselves in week 1 for “high stress” birds compared to control
Decrease in overlap to themselves in week 1 across time
Session information
## R version 4.1.1 (2021-08-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=es_ES.UTF-8
## [7] LC_PAPER=es_CR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] PhenotypeSpace_0.1.0 warbleR_1.1.27 NatureSounds_1.0.4
## [4] seewave_2.2.0 tuneR_1.4.0 ggridges_0.5.3
## [7] posterior_1.2.1 ggrepel_0.9.1 cowplot_1.1.1
## [10] tidybayes_3.0.2 modelr_0.1.8 tidyr_1.2.0
## [13] forcats_0.5.1 purrr_0.3.4 dplyr_1.0.8
## [16] lme4_1.1-29 Matrix_1.3-4 brms_2.17.0
## [19] Rcpp_1.0.8.3 GGally_2.1.2 sp_1.4-6
## [22] MASS_7.3-54 formatR_1.12 knitr_1.38
## [25] kableExtra_1.3.4 ggplot2_3.3.5 viridis_0.6.2
## [28] viridisLite_0.4.0 pbapply_1.5-0 readxl_1.4.0
## [31] lubridate_1.8.0 remotes_2.4.2
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.2 tidyselect_1.1.2 fftw_1.0-7
## [4] htmlwidgets_1.5.4 grid_4.1.1 munsell_0.5.0
## [7] codetools_0.2-18 DT_0.22 miniUI_0.1.1.1
## [10] withr_2.5.0 Brobdingnag_1.2-7 spatstat.random_2.2-0
## [13] colorspace_2.0-3 highr_0.9 rstudioapi_0.13
## [16] stats4_4.1.1 dtw_1.22-3 tensor_1.5
## [19] bayesplot_1.9.0 labeling_0.4.2 rstan_2.21.5
## [22] polyclip_1.10-0 farver_2.1.0 bridgesampling_1.1-2
## [25] coda_0.19-4 vctrs_0.4.1 generics_0.1.2
## [28] xfun_0.30 R6_2.5.1 markdown_1.1
## [31] HDInterval_0.2.2 gamm4_0.2-6 projpred_2.0.2
## [34] bitops_1.0-7 spatstat.utils_2.3-0 reshape_0.8.9
## [37] assertthat_0.2.1 promises_1.2.0.1 scales_1.2.0
## [40] gtable_0.3.0 processx_3.5.3 goftest_1.2-3
## [43] rlang_1.0.2 systemfonts_1.0.4 splines_4.1.1
## [46] spatstat.geom_2.4-0 broom_0.8.0 checkmate_2.0.0
## [49] inline_0.3.19 yaml_2.3.5 reshape2_1.4.4
## [52] abind_1.4-5 threejs_0.3.3 crosstalk_1.2.0
## [55] backports_1.4.1 httpuv_1.6.5 tensorA_0.36.2
## [58] tools_4.1.1 ellipsis_0.3.2 spatstat.core_2.4-2
## [61] raster_3.5-15 jquerylib_0.1.4 RColorBrewer_1.1-3
## [64] proxy_0.4-26 plyr_1.8.7 base64enc_0.1-3
## [67] RCurl_1.98-1.6 ps_1.6.0 prettyunits_1.1.1
## [70] rpart_4.1-15 deldir_1.0-6 zoo_1.8-10
## [73] cluster_2.1.2 magrittr_2.0.3 ggdist_3.1.1
## [76] colourpicker_1.1.1 mvtnorm_1.1-3 matrixStats_0.62.0
## [79] shinyjs_2.1.0 mime_0.12 evaluate_0.15
## [82] arrayhelpers_1.1-0 xtable_1.8-4 shinystan_2.6.0
## [85] gridExtra_2.3 rstantools_2.2.0 compiler_4.1.1
## [88] tibble_3.1.6 crayon_1.5.1 minqa_1.2.4
## [91] StanHeaders_2.21.0-7 htmltools_0.5.2 mgcv_1.8-36
## [94] later_1.3.0 RcppParallel_5.1.5 DBI_1.1.1
## [97] boot_1.3-28 permute_0.9-7 cli_3.2.0
## [100] parallel_4.1.1 igraph_1.3.0 pkgconfig_2.0.3
## [103] signal_0.7-7 terra_1.5-21 spatstat.sparse_2.1-1
## [106] xml2_1.3.3 svUnit_1.0.6 dygraphs_1.1.1.6
## [109] svglite_2.1.0 bslib_0.3.1 webshot_0.5.3
## [112] rvest_1.0.2 stringr_1.4.0 distributional_0.3.0
## [115] callr_3.7.0 digest_0.6.29 vegan_2.6-2
## [118] spatstat.data_2.2-0 rmarkdown_2.13 cellranger_1.1.0
## [121] shiny_1.7.1 gtools_3.9.2 rjson_0.2.21
## [124] nloptr_2.0.0 lifecycle_1.0.1 nlme_3.1-152
## [127] jsonlite_1.8.0 fansi_1.0.3 pillar_1.7.0
## [130] lattice_0.20-44 loo_2.5.1 fastmap_1.1.0
## [133] httr_1.4.2 pkgbuild_1.3.1 glue_1.6.2
## [136] xts_0.12.1 shinythemes_1.2.0 stringi_1.7.6
## [139] sass_0.4.1